We present a rice-wheat systen spatial decisionb support tool which relies on data from the matched landscape crop assessment surveys in eastern India and a multivariate geoadditive Bayesian model to predict entry points for system optimization in the area of interest.
2 Agronomic Functional Relationships: Crop Specific Regressions
## Also need to transform to factor for## setting up the MRF smooth.Irrig_Rev_rice_wheat$District <-as.factor(Irrig_Rev_rice_wheat$a_q103_district)## Now note that not all regions are observed,## therefore we need to remove those regions## from the penalty matrixrn <-rownames(K)lv <-levels(Irrig_Rev_rice_wheat$District)i <- rn %in% lvK <- K[i, i]set.seed(321)library(bamlss)b_rice_yield_MRF <-bamlss(f_rice_yield_MRF, data = Irrig_Rev_rice_wheat, family ="gaussian")
## First, note that we have the structured id = 'mrf1' and unstructured## spatial effect id = 're2', also indicated in the model summarysummary(b_rice_yield_MRF)
## First, note that we have the structured id = 'mrf1' and unstructured## spatial effect id = 're2', also indicated in the model summarysummary(b_wheat_yield_MRF)
# Plot the nonlinear effectplot(b_wheat_yield_MRF, model ="mu", term ="s(sowdate_fmt_wheat_day)")
Code
## Now, to predict the spatial effects we set up new data.nd <-data.frame("District"=unique(Irrig_Rev_rice_wheat$District))## Predict for the structured spatial effects.p_str_wheat_yield_MRF <-predict(b_wheat_yield_MRF, newdata = nd, term ="s(District,id='mrf1')", intercept =FALSE)## And the unstructured spatial effect.p_unstr_wheat_yield_MRF <-predict(b_wheat_yield_MRF, newdata = nd, term ="s(District,id='re2')", intercept =FALSE)## MRF smooth effect.plotmap(India_aoi_sp_bnd,x = p_str_wheat_yield_MRF$mu, id = nd$District,main =expression(mu), title ="Structured spatial effect")
Code
plotmap(India_aoi_sp_bnd,x = p_str_wheat_yield_MRF$sigma, id = nd$District,main =expression(sigma), title ="Structured spatial effect")
Code
# Plotting using bivariate maplibrary(sf)India_aoi_sf <-st_read("shp/India_aoi_sf_sp.shp")
Reading layer `India_aoi_sf_sp' from data source
`C:\Users\MMKONDIWA\OneDrive - CIMMYT\Documents\GitHub\Rice_Wheat_System_SDS_Tool\shp\India_aoi_sf_sp.shp'
using driver `ESRI Shapefile'
Simple feature collection with 47 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 82.42947 ymin: 24.28704 xmax: 88.29192 ymax: 27.85072
Geodetic CRS: +proj=longlat +datum=WGS84 +no_defs
Code
p_str_wheat_yield_MRF_dt <-as.data.frame(p_str_wheat_yield_MRF)lv_dist <-as.data.frame(lv)lv_dist_str <-cbind(lv_dist, p_str_wheat_yield_MRF_dt)India_aoi_sf_selected_str <-merge(India_aoi_sf, lv_dist_str, by.x ="NAME_2", by.y ="lv")library(biscale)library(ggplot2)data <-bi_class(India_aoi_sf_selected_str, x ="mu", y ="sigma", style ="quantile", dim =3)table(data$bi_class)
## Structured vs unstructuredIndia_aoi_sf_selected_str_dt <-subset(India_aoi_sf_selected_str, select =c("NAME_2", "mu", "sigma"))India_aoi_sf_selected_str_dt$geometry <-NULLIndia_aoi_sf_selected_unstr_str <-merge(India_aoi_sf_selected_unstr, India_aoi_sf_selected_str_dt, by ="NAME_2")library(biscale)library(ggplot2)data <-bi_class(India_aoi_sf_selected_unstr_str, x ="mu", y ="mu_unstr", style ="quantile", dim =3)table(data$bi_class)
1-1 1-2 2-1 2-2 2-3 3-1 3-2 3-3
4 6 4 2 4 2 2 6
Code
previous_theme <-theme_set(theme_bw())# create mapmap <-ggplot() +geom_sf(data = data, mapping =aes(fill = bi_class), color =NA, size =0.1, show.legend =FALSE) +bi_scale_fill(pal ="GrPink", dim =3) +labs(title ="Structured and unstructured spatial effect")bi_theme()
nd <-data.frame("District"=unique(Irrig_Rev_rice_wheat$District))# Focusing on the cross-equation correlations## Predict for the structured spatial effects.p_str_multivariate_geo_sow_MRF_corr_endo_rice_y <-predict(multivariate_geo_sow_MRF_corr_endo, newdata = nd, term ="s(District,id='mrf1')", intercept =FALSE)p_str_multivariate_geo_sow_MRF_corr_endo_rice_ydt <-as.data.frame(p_str_multivariate_geo_sow_MRF_corr_endo_rice_y)## And the unstructured spatial effect.p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y <-predict(multivariate_geo_sow_MRF_corr_endo, newdata = nd, term ="s(District,id='re2')", intercept =FALSE)# Rice sowing spatial equationplotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu1, id = nd$District, main =expression(mu(rice_sowing)), title ="Structured spatial effect")
Code
plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu1, id = nd$District, main =expression(mu(rice_sowing)), title ="Unstructured spatial effect")
Code
# Rice yield spatial equationplotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu2, id = nd$District, main =expression(mu(rice_yield)), title ="Structured spatial effect")
Code
plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu2, id = nd$District, main =expression(mu(rice_yield)), title ="Unstructured spatial effect")
Code
# Wheat sowing equationplotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu3, id = nd$District, main =expression(mu(wheat_sowing)), title ="Structured spatial effect")
Code
plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu3, id = nd$District, main =expression(mu(wheat_sowing)), title ="Unstructured spatial effect")
Code
# Wheat yield spatial equationplotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu4, id = nd$District, main =expression(mu(wheat_yield)), title ="Structured spatial effect")
Code
plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu4, id = nd$District, main =expression(mu(wheat_yield)), title ="Unstructured spatial effect")
Code
# Focusing on the cross-equation correlations## MRF smooth effect.plotmap(India_aoi_sp_bnd,x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$lambda24, id = nd$District,main =expression(lambda(rice, wheat)), title ="Structured spatial effect")
Code
## Random effects.plotmap(India_aoi_sp_bnd,x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$lambda24, id = nd$District, main =expression(lambda(rice, wheat)), title ="Unstructured spatial effect")
Code
# Rice sowing equation : Non-linear relationships# s(gw_2017) + s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + s(sd_onset_82_15)plot(multivariate_geo_sow_MRF_corr_endo, model ="mu1", term ="s(gw_2017)")
Code
plot(multivariate_geo_sow_MRF_corr_endo, model ="mu1", term ="s(onset_2017)")
Code
plot(multivariate_geo_sow_MRF_corr_endo, model ="mu1", term ="s(monsoon_onset_dev)")
Code
plot(multivariate_geo_sow_MRF_corr_endo, model ="mu1", term ="s(median_onset_82_15)")
Code
plot(multivariate_geo_sow_MRF_corr_endo, model ="mu1", term ="s(sd_onset_82_15)")
Code
# Rice yield equation# s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice)plot(multivariate_geo_sow_MRF_corr_endo, model ="mu2", term ="s(Res_rice_sow)")
Code
plot(multivariate_geo_sow_MRF_corr_endo, model ="mu2", term ="s(g_q5305_irrig_times_rice)")
Code
plot(multivariate_geo_sow_MRF_corr_endo, model ="mu2", term ="s(nperha_rice)")
Code
plot(multivariate_geo_sow_MRF_corr_endo, model ="mu2", term ="s(p2o5perha_rice)")
Code
# Wheat sowing equation# s(harvest_day_rice) + s(gw_2018)plot(multivariate_geo_sow_MRF_corr_endo, model ="mu3", term ="s(harvest_day_rice)")
Code
# (multivariate_geo_sow_MRF_corr_endo, model = "mu3", term = "s(gw_2018)")# Wheat yield equationplot(multivariate_geo_sow_MRF_corr_endo, model ="mu4", term ="s(Res_wheat_sow)")
Code
# Fitted valuesmultivariate_geo_sow_MRF_corr_endo_fitted_values <- multivariate_geo_sow_MRF_corr_endo$fittedmultivariate_geo_sow_MRF_corr_endo_fitted_values <-as.data.frame(multivariate_geo_sow_MRF_corr_endo_fitted_values)# Merge the fitted results to the data and exportIrrig_Rev_rice_wheat_Mult_Res <-cbind(Irrig_Rev_rice_wheat, multivariate_geo_sow_MRF_corr_endo_fitted_values)summary(lm(sowdate_fmt_rice_day ~ mu1, Irrig_Rev_rice_wheat_Mult_Res))
Call:
lm(formula = sowdate_fmt_rice_day ~ mu1, data = Irrig_Rev_rice_wheat_Mult_Res)
Residuals:
Min 1Q Median 3Q Max
-44.574 -5.397 -0.269 4.986 44.985
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.03383 3.41385 -0.01 0.992
mu1 1.00018 0.01771 56.48 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.093 on 4525 degrees of freedom
Multiple R-squared: 0.4135, Adjusted R-squared: 0.4134
F-statistic: 3190 on 1 and 4525 DF, p-value: < 2.2e-16
## Also need to transform to factor for## setting up the MRF smooth.Irrig_Rev_rice_wheat$District <-as.factor(Irrig_Rev_rice_wheat$a_q103_district)## Now note that not all regions are observed,## therefore we need to remove those regions## from the penalty matrixrn <-rownames(K)lv <-levels(Irrig_Rev_rice_wheat$District)i <- rn %in% lvK <- K[i, i]# First stage analyticslibrary(sp)library(mgcv)library(bamlss)# Multivariate geoadditive model# remotes::install_git("https://git.uibk.ac.at/c4031039/mvnchol")# install_github("https://github.com/meteosimon/mvnchol")library(mvnchol)# Rice first stagef_sow_rice_1st_stage <-list( sowdate_fmt_rice_day ~1+ rice_duration_class_long +s(gw_2017) +s(May_tmax_17) +s(onset_2017) +s(monsoon_onset_dev) +s(median_onset_82_15) +s(sd_onset_82_15) +s(District, bs ="mrf", xt =list("penalty"= K)) +s(District, bs ="re"))f_sow_rice_1st_stage_MRF <-bamlss(f_sow_rice_1st_stage, data = Irrig_Rev_rice_wheat, family ="gaussian")